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We present detailed molecular dynamics results for the displacement autocor- 
relation spectra of the Peyrard-Bishop model of thermal DNA denaturation. As 
the phase transition is approached, the spectra depend on whether the wavelength 
is smaller than, or exceeds the correlation length. In the first case, the spectra 
are dominated by a single peak, whose frequency approaches the bare acoustic 
frequency of the harmonic chain, and whose linewidth approaches zero as T c — T. 
In the second case, a central peak (CP) feature is dominant, accounting for most 
of the weight; the linewidth of the CP appears to be temperature- independent. 
We also present force autocorrelation spectra which may be relevant for analyzing 
the statistical properties of localized modes. 

1 Introduction 

The thermal denaturation of DNA, i.e. the separation of the two strands 
upon heating, is a typical thermodynamic instability. It can be modelled 
along the lines of other thermodynamic instabilities (e.g. wetting, solid-on- 
solid adsorption), by associating a single, one-dimensional coordinate with 
the distance of a base pair[Q. Details can be found in this volume [Q] and 
in the original literature cited there. 
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The equilibrium properties of the system near the phase transition are 
characterized by a divergent correlation length £, and a discontinuity in the 
specific heat; in other words, this is a second-order transition; the feature 
which sets it apart from other structural, or order-disorder transitions is 
that, as the transition temperature is approached from below, the order 
parameter diverges, i.e. the low-temperature phase becomes continuously 
unstable. 

In this paper, we present results for the dynamical correlations of the 
order parameter, obtained by numerical simulation. At low and interme- 
diate temperatures, the spectra appear to be dominated by the properties 
of localized anharmonic motion ( "discrete breathers" ) . As the critical tem- 
perature is approached from below, the spectra depend solely on whether 
the wavelength is smaller or larger than the correlation length. In the first 
case, they reflect the dynamics of "islands" of the high temperature phase. 
In the second case, they are dominated by a strong central peak, whose 
width appears to depend on the wavevector, but not on the temperature. 

The paper is structured as follows: Section 2 introduces the notation and 
numerical procedure. Section 3 presents the main results, and an analysis 
along the lines of relaxational/oscillational phenomenology. Section 4 is a 
sketch of a tentative, alternative theory, along the lines of the Mori-Zwanzig 
projection operator formalism. Section 5 presents a brief summary and 
discussion. 



2 Notation and numerical procedure 

We consider the "minimal" Hamiltonian model of homogeneous DNA de- 
naturation proposed by Peyrard and Bishop]]]] (PB), 



H(y) = J2 



v 1 

■Y + ^Vn-Vn-lf +V{y n ) 



(1) 



where y n ,Pn are dimensionless, canonically conjugate coordinates and mo- 
menta of the nth base pair transverse to the chain, and V(y) = (1 — e~ v ) 2 . 
R is a dimensionless parameter which describes the relative strength of 
on-site vs. elastic interactions; here R = 10.1. 

The thermodynamic properties of (|l|) have been reviewed in Ref. 2. This 
work describes the spectra of dynamical correlations 

S AA (q,iu) =- — e --*£V<M«-™) < A n (t)A rn (0) > (2) 

00 m.n 

where, in this paper, mostly A n = y n . The integral of @ over all frequen- 
cies, 

/OO 1 
dio S AA {q, u) = - £ e lqa{n - m) < A n A m > (3) 
-OO ™ 
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can be computed exactly using the transfer-integral result for the equal time 
correlations (cf. Eq. 55 of Ref. 2). It is expedient to consider normalized 
spectral functions 

The angular brackets in Eqs. (||)-(||) denote canonical ensemble averages. 
Typically, we implemented this by repeated molecular dynamics (MD) 
simulations of the system for many different initial conditions, Fourier- 
transforming the spatiotemporal correlations obtained from each run, and 
averaging over all runs to obtain the final result. 
The equations of motion, 

Vn = (Vn+1 + Vn-1 ~ tyn) - V (y n ) ; 11=1,2, ...N (5) 

with periodic boundary conditions, ya = j/jv, Un+i = yi, and typical sys- 
tem size N — 1024, were numerically integrated for an interval T = 410, 
using a 4-th order Runge-Kutta algorithm, with a time step equal to 0.02. 
Initial conditions were "canonical" , in the sense that (i) the velocities y n 
were Gaussian-distributed, and (ii) the positions y n were random variables 
distributed according to the potential energy part of the Hamiltonian (|l|) ; 
in addition, the system was "thermalized" for a certain time, using a Nose 
procedure^). 



3 Spectra: Phonons vs. central peak 

At intermediate temperatures, the spectra are characterized by an an- 
harmonic phonon component and a strong, low-frequency intensity (cf. 
Fig. 0a); this low- frequency component becomes even more pronounced at 
lower values of the wavevector. There is considerable residual structure in 
the spectrum; in particular, a secondary peak at lower frequencies appears 
to be a consistent feature. 

As the temperature increases, and the instability approaches, the struc- 
ture becomes significantly simpler. The decisive quantity is the correlation 
length £. If the wavelength is shorter than £, i.e. q£/ (2ir) > 1, the spectrum 
in effect probes the "droplets" of the high-temperature phase, of typical 
size £, which are present in the low-temperature phase; consequently, the 
main feature of Fig. [l]b is a peak, from the acoustic phonons Q. At the 
smallest values of the ratio g£/(27r), a central peak (CP) feature begins to 
grow, and eventually dominates the spectrum at values g£/(27r) << 1; this 
is the case in Fig. ||. 

A first attempt to analyze the data can be made in terms of a phe- 
nomenological relaxation/oscillation spectral function, similar to the one 
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Figure 1: Normalized dynamical correlation spectra. Panel (a) presents 
the results at T — 0.5, for selected values of q. The inset shows the zero- 
frequency details (linear frequency scale). Panel (b) presents results at 
T = 1.05, for a variety of q values. The inset shows the ratio of the cor- 
relation length £ and the wavelength. If the wavelength is smaller than 
the correlation length, the spectrum is dominated by the phonon peak. At 
this length scale, the spectra probe the "droplets" of the high-temperature 
phase present in the system. Note the gradual buildup of a central peak at 
the lowest values of q. 
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Figure 2: Normalized dynamical correlation spectra at T — 0.7, and q/n = 
0.07. At long-wavelengths (compared to £) the spectrum is dominated by 
the CP feature. The inset shows that the tail of the spectra drops off with 
a different slope. The fit has been obtained using Eq. @. 



used in analyzing structural phase transitions) 
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where T = Tq + S 2 /(j — iu>) is a relaxational memory kernel, and the in- 
dependence of all the parameters has been suppressed. If To <C 5 2 /'~f, it is 
possible for the spectrum (^) to split into phonon-like 

^(^).l(l-p) ( ^ 2 _^ (cjro)2 , .» 7 (7) 

and CP 

i y 

S yy (q^)^- P ^ 2+ ^ 2 , C^«7 (8) 

contributions, where uj 2 ^ — uj 2 + S 2 , p — 5 2 /oj 2 ^, and 7' = 7(1 — p). 

Fig. ||a shows that the CP linewidth 7' is largely independent of tem- 
perature; it does however depend on q, roughly linearly, as long as q^/(2ir) 
does not exceed unity, i.e. as long as the CP is appreciable. 



S yy (q,uj) = Im— — (6) 

ituj <j)i{ — v 1 — lUJl 
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Figure 3: Linewidths. Panel (a) presents the linewidth of the central peak 
as a function of wavelength. The dashed line has unit slope. Panel (b) 
presents the linewidth of the q = ir phonon as a function of temperature. 
The linewidth extrapolates to zero at a temperature not far from T c . 
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Returning to the phonon peak, it is possible to follow the decreasing 
phonon linewidth as the transition is approached and the dynamics - with 
the exception of the very long wavelengths g£/(27r) < 1 - evolves towards 
the harmonic limit(Fig. |3|b). Our data is consistent with a linear critical 
slowing down, i.e. T p h oc T c — T . 

4 Force autocorrelations 

It is instructive to consider the spectrum of the force autocorrelations, i.e. 
A n = f n = —V'(y n ) in Eqs. - (||). At low temperatures one might try to 
describe the spectra in terms of a crude model of independent local modes 
(ILM), i.e. site-independent solutions of the R — > oo (anticontinuum) limit 
of (ph. The explicit form of the second time derivative is 



where < A < X max — V2, < S < 2tt depend on the initial conditions 
and e\ = 1 — A 2 /2 is the energy of the ILM. Since R ^> 1, the above form 
should be a good approximation to exact one-site discrete breathers Q. The 
canonical average of autocorrelations of (0) is 



where Q(A) = Z^ 1 exp (—e\/T) and Z is determined from the normaliza- 
tion condition J dX Q(X) = 1. The MD force spectra at T = 0.2 are 
shown in Fig. ^a, along with a numerical Fourier transform of (|l0|). Over- 
all agreement is satisfactory, except for the very low frequency part of the 
spectrum. 

At higher temperatures, force autocorrelation spectra exhibit the fol- 
lowing features: (i) a very pronounced dip occurring almost exactly at 
the bare acoustic frequencies uj q — (2/%/i?) sin(q/2), characteristic of the 
high-temperature phase, and (ii) a roughly ^-independent decay at higher 
frequencies. The observed form of the force spectra motivates an improved 
version of (||), with = T/S(q), as demanded by the Mori-Zwanzig 
projection operator formalism]^], and ReT = Fo exp(— lo/uj c ) + Ti, where 
Pi = a(l — uj/tbqf/ 2 if uj < Lu q , and Ti — b(uj~ ojq) 1 ^ 2 exp(— lu/uj c ) if lu > uj q . 
Preliminary fits obtained with this improved Ansatz for the memory func- 
tion (approximating ImT by a constant, equal to its value at the peak, and 
setting lu c equal to the value obtained by the exponental decay constant of 
force spectra, cf. above and Fig. ||b) are shown in Fig. [| they seem to 
reproduce the MD data much better, using the same number of adjustable 
parameters. 



fx,s(t) = A 2 e 



2 1/2 cos(Xt + 6)- e A ' 



(9) 



1 - e\ /2 cos(A£ + 5) 




(10) 
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Figure 4: Normalized spectra of force autocorrelations. Panel (a) shows the 
results at low temperatures. The wavy curve is a theoretical estimate based 
on an independent localized mode picture, i.e. the normalized spectrum 
of (jio|). Panel (b) presents results at higher temperatures. Note (i) the 
extreme dip, almost to zero intensity, which occurs almost exactly at the 
acoustic phonon frequencies (inset), and (ii) the exponential decay at higher 
frequencies. 
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Figure 5: The spectra S yy at T = 1.0. Fits are obtained with the im- 
proved Mori-Zwanzig Ansatz, which incorporates a memory kernel with a 
dip, similar to the one shown in Fig. ^b. 

5 Concluding remarks 

The MD data presented show that the critical dynamics of the Peyrard- 
Bishop model of DNA thermal denaturation can be thought of as follows: 
At length scales shorter than the correlation length £, which correspond to 
"droplets" of the high temperature phase, the system reflects the properties 
of the unstable phase; oscillatory dynamics of the soft, acoustic phonons 
is the result. The linewidth of these phonons appears to vanish linearly as 
T c — T{ "critical slowing down" ) . At length scales longer than the correla- 
tion length, the dynamics is dominated by the central peak. Fluctuations 
are stronger, as evidenced from the divergence of the static structure fac- 
tor S(q); the typical time scales of these fluctuations appears however to 
be non-critical, ft remains a challenge to the theory to establish whether 
these "non-critical" , g-dependent dynamics can be associated with localized 
excitations. The preliminary analysis performed at low temperatures sug- 
gests that a picture of independent localized modes provides a reasonable 
description of the force autocorrelations - with the important exception of 
the very low frequency regime. Perhaps a more detailed theory of discrete 
breather statistical mechanics can improve our understanding of this part 
of the spectra. 
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